Trimmed and filtered reads were aligned using STAR (v2.7.0f). Read counting was performed using featureCounts (v1.6.2) with parameters “-s 1 -Q 255”. Gene annotation was generated by merging GENCODE (vM21), 3’UTR annotations in RefSeq (mm10, UCSC), and de novo 3’UTR sites discovered from the data.
Count table is read into R. A sample condition table is generated by parsing filenames.
# gene count file
counts_file <- "datasets/pit_utrseq_gencode_vM21_refSeq201906.txt"
# Obtain gene count table and sample conditions
temp <- get_count_table(counts_file, fixname = T)
count_table <- temp[[1]]
sampleCondition <- temp[[2]]
gene_info <- temp[[3]]
rm(temp)
colnames(count_table) <- gsub("p", "P", colnames(count_table))
# Get all chromosome Y genes and 3 chrX genes that are consistently female-biased
chrXYgenes <- as.character(gene_info[grepl("chrY", gene_info[,2]),1])
chrXYgenes <- c(chrXYgenes, "ENSMUSG00000086503.3", "ENSMUSG00000086503.4", "ENSMUSG00000099312.1")
chrMgenes <- as.character(row.names(gene_info[grepl("chrM", gene_info[,2]),]))
We used RUVseq to remove unwanted variables with empirically identified negative control genes.
# constructc edger expression set to obtain cpm values
y <- DGEList(counts=count_table, genes=row.names(count_table), group=sampleCondition$group)
y <- calcNormFactors(y)
count_cpm <- cpm(y)
# get spike-in gene names
spikes <- rownames(count_table)[grep("^ERCC", row.names(count_table))]
# filter for only genes expressed at cpm > 2 in at least 10 samples. Also removed mitochondrial genes and spike-ins
count_table_fil <- count_table[which(rowSums(count_cpm > 2) >= 10), ]
count_table_fil <- count_table_fil[!row.names(count_table_fil) %in% chrMgenes,]
count_table_fil <- count_table_fil[!row.names(count_table_fil) %in% spikes,]
# construct edger expression set again with filtered data
y_fil <- DGEList(counts=count_table_fil, genes=row.names(count_table_fil), group=sampleCondition$age)
y_fil <- calcNormFactors(y_fil)
# get log CPM counts
logCPM <- cpm(y_fil, log=TRUE, prior.count=1)
# Correct batch effect/unwanted variables using RUV seq
row.names(sampleCondition) <- colnames(count_table_fil)
set <- newSeqExpressionSet(as.matrix(count_table_fil),
phenoData = sampleCondition)
# upper quartile normalzation
set <- betweenLaneNormalization(set, which="upper")
# identify empirical control genes by performing an DE analysis across all conditions
design <- model.matrix(~group, data=pData(set))
y_fil <- DGEList(counts=count_table_fil, genes=row.names(count_table_fil), group=sampleCondition$group)
y_fil <- calcNormFactors(y_fil)
y_fil <- estimateDisp(y_fil, robust = TRUE, design = design)
fit <- glmQLFit(y_fil, robust = T, design)
lrt <- glmQLFTest(fit, coef=2:10)
top <- topTags(lrt, n=nrow(count_table_fil))$table
top$genename <- genename[top$genes]
# identified any genes with an FDR > 0.1 as negative controls
empirical <- subset(top, FDR > 0.1)$genes
# RUVseq normliazation
set1 <- RUVg(set, empirical, k=2)
# obtain normalized counts
logCPMc <- log2(normCounts(set1)+1)
qPCR of 183 puberty-related and control genes have been performed previously using the same samples. Here, we compared the expression quantification using 3’UTRseq (log2 normalized read counts) and that using qPCR (delta Ct)
# get qPRC results
all_data_normed_info_noC <- readRDS("datasets/fluidigm_all_data_normed_info_noC.rds")
# fix gene names
all_data_normed_info_noC$Primer[all_data_normed_info_noC$Primer == "Rab7l1"] <- "Rab29"
fluidigm_genes <- unique(all_data_normed_info_noC$Primer)
all_genes <- fluidigm_genes
# reformat pcr results table
PCR_puberty <- subset(all_data_normed_info_noC, tissue == "P") %>%
mutate(Sample = gsub("_", "", Sample)) %>%
group_by(Sample, Primer, age, gender) %>%
summarise(Ct= mean(Ct)) %>%
as.data.frame()
# compare pcr results and utrseq results (before and after normalization, only post-normalization data shown)
exp_compare_cor_values <- compare_qPCR(logCPM, logCPMc, PCR_puberty, outname, method = "spearman", ifgenplot = F)
compare_sum <- compare_qPCR_summaryplot(exp_compare_cor_values$cor_values)
exp_compare_cor_values$plots_after_correction$Pd37M4
Figure 1: Figure S1C: summary of UTRseq vs. qPCR quantification comparison
compare_sum[[2]]
Figure 2: Figure 1C: example scatter plot comparing UTRseq and qPCR
Pearson correlation coefficient between normalized UTRseq counts and qPCR delta Ct ranges from 0.7 to 0.78.
PCA plot and correlation heatmap using RUVseq normalized counts (log2 transformed).
# get PCA plots
logCPMc_PCAplot <- gen_pca_plot(logCPMc, npc=20, sample_label=F, label =F)
logCPMc_PCAplot[[1]]
Figure 3: Figure 1D: PCA plot of all samples
hmcol <- colorRampPalette(c("orange", "white", "darkmagenta"))(250)
plotHM(logCPMc, cols = hmcol, new_col_names = sapply(strsplit(colnames(logCPMc), "_"), "[[", 2), show_col_names = F, clust_method = "complete", anno = sampleCondition)
Figure 4: Figure S2A: correlation heatmap of samples
# calculate mean person correlation between biological reps
sample_cors <- cor(logCPMc)
sample_cons <- unique(sampleCondition$group)
x <- sample_cons[1]
sample_cors_means <- lapply(sample_cons, function(x) {
cors <- sample_cors[grepl(x, row.names(sample_cors)), grepl(x, colnames(sample_cors))]
cors[cors==1] <- NA
return(list(mean(cors, na.rm=T), cors))}
)
mean_cor <- round(sapply(sample_cors_means,"[[", 1),2)
Average pearson correlation between biological replicates (samples from the same biological condition) ranges from 0.95 to 0.97
Pairwise differential analysis was performed using EdgeR. Genes with a FDR < 0.05 and an absolute fold change> 1.5 are considered significant.
Comparisons are performed between male and female samples at each age.
all_detected_genes <- genename[row.names(count_table_fil)]
y_fil <- DGEList(counts=count_table_fil, genes=row.names(count_table_fil), group=sampleCondition$group)
y_fil <- calcNormFactors(y_fil)
design <- model.matrix(~0+group+W_1+W_2, pData(set1))
y_fil <- estimateDisp(y_fil, robust = TRUE, design = design)
fit <- glmQLFit(y_fil, design=design, robust = TRUE)
# compare between sexes at each age.
pair_contrasts <- makeContrasts(d12_sex = groupd12F - groupd12M,
d22_sex = groupd22F - groupd22M,
d27_sex = groupd27F - groupd27M,
d32_sex = groupd32F - groupd32M,
d37_sex = groupd37F - groupd37M,
levels = design)
get_test_result <- function(contrasts=pair_contrasts, name){
qlf <- glmQLFTest(fit, contrast = contrasts[, name])
topTg <- topTags(qlf, n=nrow(y_fil$counts))
de <- as.data.frame(topTg[[1]])
de$genename <- genename[de$genes]
return(de)
}
de_result_list <- lapply(colnames(pair_contrasts), function(x) get_test_result(pair_contrasts, x))
names(de_result_list) <- colnames(pair_contrasts)
de_sig_list <- lapply(de_result_list, function(x) subset(x, FDR < 0.05 & abs(logFC) > log2(1.5)))
names(de_sig_list) <- colnames(pair_contrasts)
templist <- flatten(lapply(de_sig_list, function(x) list(subset(x, logFC>0), subset(x, logFC<0))))
names(templist) <- flatten(lapply(names(de_sig_list), function(x){
sapply(c("up", "down"), function(y) paste(x, y, sep = "_"))
}))
sig_genes_names_list_direction <- lapply(templist, "[[", "genename")
sig_genes_names_list <- lapply(de_sig_list, "[[", "genename")
sex_sig_genes_df <- get_sig_genes_df_from_list(de_sig_list)
colnames(sex_sig_genes_df) <- gsub("_sex_up", "_F", colnames(sex_sig_genes_df))
colnames(sex_sig_genes_df) <- gsub("_sex_down", "_M", colnames(sex_sig_genes_df))
upset(sex_sig_genes_df, sets=rev(c("d12_F", "d22_F", "d27_F", "d32_F","d37_F", "d12_M", "d22_M", "d27_M","d32_M", "d37_M")), keep.order = T,
order.by = "freq",
queries = list(list(query = intersects, params = list("d27_F", "d32_F", "d37_F"), color="tomato", active=T),
list(query = intersects, params = list("d27_M","d32_M", "d37_M"), color="steelblue", active=T)))
Figure 5: Figure 4A: overlap between sex-biased genes
# get sex biased genes in 2 ages
temp <- table(unlist(sig_genes_names_list_direction[c("d27_sex_up", "d32_sex_up","d37_sex_up")]))
f_biased_genes_used <- names(temp[temp !=1])
f_biased_genes_3 <- names(temp[temp ==3])
temp <- table(unlist(sig_genes_names_list_direction[c("d27_sex_down", "d32_sex_down","d37_sex_down")]))
m_biased_genes_used <- names(temp[temp !=1])
m_biased_genes_3 <- names(temp[temp ==3])
# f_biased_genes_used and m_biased_genes_used are used as input for Lisa
# sex age interaction
sex_age_contrasts <- makeContrasts(d12_d22 = (groupd22F - groupd12F) - (groupd22M - groupd12M),
d22_d27 = (groupd27F - groupd22F) - (groupd27M - groupd22M),
d27_d32 = (groupd32F - groupd27F) - (groupd32M - groupd27M),
d37_d32 = (groupd37F - groupd32F) - (groupd37M - groupd32M),
d32_d22 = (groupd32F - groupd22F) - (groupd32M - groupd22M),
d37_d27 = (groupd37F - groupd27F) - (groupd37M - groupd27M),
d37_d22 = (groupd37F - groupd22F) - (groupd37M - groupd22M),
levels = design)
de_sex_age_result_list <- lapply(colnames(sex_age_contrasts), function(x) get_test_result(sex_age_contrasts, x))
names(de_sex_age_result_list) <- colnames(sex_age_contrasts)
de_sex_age_sig_list <- lapply(de_sex_age_result_list, function(x) subset(x, FDR < 0.05 & abs(logFC) > log2(1.5)))
sig_sex_age_genes_names_list <- lapply(de_sex_age_sig_list, function(x) x$genename)
#sig_genes_puberty <- lapply(sig_genes_names_list, function(x) x[x %in% all_puberty_genes])
d12_d22_sex_diff <- de_sex_age_sig_list$d12_d22$genename
#factor_plot(logCPMc, sort(d12_d22_sex_diff), ncols=5, name=F)
sig_df <- melt(sig_genes_names_list) %>%
separate(L1,into = c("age",NA))
# plot selected d12-d22 sex biased genes
ave_expr <- rowMeans(logCPMc)
names(ave_expr) <- genename[names(ave_expr)]
# figure2 genes are selected based on 1) d12 sex biased genes , 2) gene function mentioned in the text and 3) expression levels
figure2_genes1 <- c("Chrna4","Kcna4","Th","Drd4", "Asb4", "Fshb","Lhb","Serpine2", "Nupr1","Gpr101")
factor_plot(logCPMc, figure2_genes1, ncols = 5, sig_df = sig_df)
# get top m-/f-biased genes that are consistent from PD27
top_fbiased <- sort(ave_expr[names(ave_expr) %in% f_biased_genes_3], decreasing = T)[1:4]
top_mbiased <- sort(ave_expr[names(ave_expr) %in% m_biased_genes_3], decreasing = T)[1:4]
factor_plot(logCPMc, c(names(top_fbiased), names(top_mbiased)), ncols = 5, sig_df = sig_df)
enrich_list_direction <- lapply(sig_genes_names_list_direction, function(x) gprofiler(x, organism = "mmusculus", custom_bg = all_detected_genes, hier_filtering = "moderate", src_filter = c("GO", "KEGG","REAC","CORUM","HP","HPA","OMIM"),min_set_size = 5, max_set_size = 5000,min_isect_size = 3))
male_biased_pathways <- enrich_list_direction[c("d27_sex_down", "d32_sex_down", "d37_sex_down")]
female_biased_pathways <- enrich_list_direction[c("d27_sex_up", "d32_sex_up", "d37_sex_up")]
maleenrichmentmap <- enrichmap_like_pie(male_biased_pathways, layout_used = "kk", radius_scale = 0.05)
femaleenrichmentmap <- enrichmap_like_pie(female_biased_pathways, layout_used = "kk", radius_scale = 0.07)
maleenrichmentmap[[2]]
Figure 6: Figure 2F: enrichment map-like representation of pathway enrichment of male genes at PD27, 32, and 37
maleenrichmentmap[[3]]
Figure 7: Figure 2F: enrichment map-like representation of pathway enrichment of male genes at PD27, 32, and 37
femaleenrichmentmap[[2]]
Figure 8: Figure 2G: enrichment map-like representation of pathway enrichment of female genes at PD27, 32, and 37
femaleenrichmentmap[[3]]
Figure 9: Figure 2G: enrichment map-like representation of pathway enrichment of female genes at PD27, 32, and 37
The genes shown in this heatmap includes:
that are not already shown in Figure 2
# get sex biased genes in 2/3 ages
sex_biased_genes_to_show <- unique(c(unlist(sig_genes_names_list[1:2]),d12_d22_sex_diff))
# exclude genes that are shown in figure 2 already
sex_biased_genes_to_show <- sex_biased_genes_to_show[!sex_biased_genes_to_show %in% c(figure2_genes1, names(top_fbiased), names(top_mbiased))]
#sex_biased_genes_to_show <- c(sex_chr_genes, "Kcna4","Th", d22_rest)
sex_genes_sum <- table(unlist(sig_genes_names_list[1:5]))
temp_df <- data.frame(value = sex_biased_genes_to_show, stringsAsFactors = F)
sex_gene_df <- left_join(temp_df, melt(sig_genes_names_list_direction)) %>%
# filter(value %in% sex_biased_genes_to_show) %>%
separate(L1, "_", into = c("age", NA, "dir")) %>%
spread(age, dir) %>%
mutate(n = sex_genes_sum[as.character(value)]) %>%
dplyr::select(value, d12, d22, d27, d32, d37) %>%
arrange(d12, d22, d27, d32, d37) %>%
column_to_rownames("value")
#sex_gene_df$d12_d22 <- ifelse(row.names(sex_gene_df) %in% d12_d22_sex_diff, "yes","no")
anno_cols_sex <-rep(list(c("up" = "tomato","down" = "steelblue")), 5)
names(anno_cols_sex) <- c("d12","d22","d27","d32","d37")
#anno_cols_inter <- list("d12d22" = c("yes" = "darkmagenta"))
p <- plot_gene_hm(row.names(sex_gene_df),
logCPMc,
clust_cols = F,
clust_rows = F,
row_cut = 1,
anno_colors = c(anno_cols, anno_cols_sex),
anno_rows = sex_gene_df[, c("d12","d22","d27","d32","d37")],
colnames = F,
# fontsize_row = 8,
colors = hmcol_yp,
sampleCon = sampleCondition)
plot(p$gtable)
Figure 10: Figure S3A: heatmap of selected sex-biased genes
These genes are TFs that are sex-biased at 2 out of 3 ages (PD27, 32, and 37) and are not shown in Figure 3.
# plots sex-biased tfs
sex_tfs <- c(m_biased_genes_used[m_biased_genes_used %in% alltfs], f_biased_genes_used[f_biased_genes_used %in% alltfs])
temp_df <- data.frame(value = sort(sex_tfs), stringsAsFactors = F)
sex_gene_df <- left_join(temp_df, melt(sig_genes_names_list_direction)) %>%
separate(L1, "_", into = c("age", NA, "dir")) %>%
filter(!is.na(age)) %>%
# spread(age, dir) %>%
pivot_wider(names_from = "age", values_from = "dir") %>%
mutate(sex = ifelse(d27=="down"|d32=="down"|d37=="down","down","up")) %>%
group_by(sex) %>%
arrange(d27, d32, d37,.by_group = T) %>%
column_to_rownames("value")
anno_cols_sex <-rep(list(c("up" = "tomato","down" = "steelblue")), 5)
names(anno_cols_sex) <- c("d12","d22","d27","d32","d37")
p_hm_tf <- plot_gene_hm(row.names(sex_gene_df),
logCPMc,
clust_cols = F,
clust_rows = F,
row_cut = 1,
anno_colors = c(anno_cols, anno_cols_sex),
anno_rows = sex_gene_df[, c("d27","d32","d37")],
colnames = F,
# fontsize_row = 8,
colors = hmcol_yp,
sampleCon = sampleCondition)
plot(p_hm_tf$gtable)
Figure 11: Figure 3B: Expression of sex-biased TFs
# run CEMtools with default parameters
cem <- cemitool(as.data.frame(logCPMc), merge_similar = T, network_type = "signed")
# get module assignment
modules <- module_genes(cem)
modules$genename <- genename[modules$genes]
cemi_genes <- modules$genes
# transform into a list format
module_list <- lapply(split(modules, modules$modules), "[[", "genes")
module_list_genename <- lapply(split(modules, modules$modules), "[[", "genename")
# get average scaled expression for each gene for plotting in the next section
ave <- ave_count_table(logCPMc)
# get top hub genes
hubs <- CEMiTool::get_hubs(cem, n=20)
hubs_genename <- lapply(hubs, function(x) genename[names(x)])
-Left:Heatmap showing the scaled normalized expression levels of each gene, clustered separately within modules.
-Right: plot summarizing scaled normalized expression patterns of each gene in each module. Thick line represent the average profile of the module.
-Bottom: top 10 hub genes for each module
# scale normalized count data (log transformed) and reorder columns
plotdata <- t(scale(t(logCPMc)))
colnames(plotdata) <- sapply(strsplit(colnames(plotdata), "_"), "[[", 2)
plotdata <- plotdata[, c(colnames(plotdata)[grepl("M", colnames(plotdata))], colnames(plotdata)[grepl("F", colnames(plotdata))])]
get_clust_gene_order <- function(clust_num){
# function to cluster genes within each cluster
used_genes <- module_list[[clust_num]]
used_data <- plotdata[row.names(plotdata) %in% used_genes, ]
phm <- pheatmap(used_data,
clust_cols = F,
clustering_distance_rows = "correlation",
clustering_method = "ward.D2",
silent = T)
order <- used_genes[phm$tree_row$order]
return(order)
}
# reorder genes based on clustering results within each cluster
reorder <- unlist(lapply(1:9, get_clust_gene_order))
cluster_length <- sapply(module_list, length)
# column annotation data
col_anno <- pData(set1)
row.names(col_anno) <-sapply(strsplit(row.names(col_anno), "_"), "[[", 2)
# set cutoffs for extreame values
all_values <- as.vector(as.matrix(plotdata))
custome_breaks <- unique(c(min(all_values), seq(quantile(all_values, 0.005), 0, length = 125), seq(0, quantile(all_values, 0.995), length=125), max(all_values)))
# generate heatmap
phm_modulegenes <- pheatmap(plotdata[reorder, ],
cluster_rows = F,
cluster_cols = F,
show_rownames = F,
annotation_col = col_anno[, c("age", "sex")],
annotation_colors = anno_cols,
color = hmcol,
breaks = custome_breaks,
gaps_row = cumsum(cluster_length),
silent = T)
# get a summarized plot of module gene expression
ave <- ave_count_table(logCPMc)
module_assignment <- modules$modules
names(module_assignment) <- modules$genes
pclusters <- plot_clusters(module_assignment, ave, "", ncols=1) + theme(axis.text.x = element_blank(), axis.text.y = element_text(size=14), strip.text = element_text(size=14))
hub_df <- as.data.frame(sapply(hubs_genename, function(x) paste(x[1:10], collapse = ",")))
colnames(hub_df) <- "Top10_hub_genes"
# put together
p <- as_ggplot(phm_modulegenes$gtable) + pclusters +
plot_layout(widths = c(4,1.5,1)) +
theme_bw()
p
Figure 12: Figure 4:co-expression module expression and hub genes
plot(gridExtra::tableGrob(hub_df))
Figure 13: Figure 4:co-expression module expression and hub genes
top5_hub_genes <- unname(unlist(lapply(hubs_genename, function(x) x[1:5])))
factor_plot(logCPMc, top5_hub_genes, ncols = 5)
Figure 14: Figure S4: Expression patterns of top 5 hub genes
hubs_all <- bind_rows(lapply(CEMiTool::get_hubs(cem, n="all"),enframe, name = "genes", value = "gene_connectivity"), .id="modules")
modules_out <- left_join(modules,hubs_all) %>%
arrange(modules, desc(gene_connectivity)) %>%
select(modules, everything())
#write.xlsx(modules_out, "7_supp_table_co_expression_module_assignment.xlsx")
module_ME <- mod_summary(cem, "eigengene")
module_ME_df <- melt(module_ME) %>%
separate(variable, "_", into=c("id", "condition")) %>%
mutate(age=substr(condition, 2,4), sex=substr(condition, 5,5), rep=substr(condition, 6,6))
row.names(module_ME) <- module_ME$modules
module_ME <- module_ME[, -1]
p <- pheatmap(cor(t(module_ME)), silent = T)
plot(p$gtable)
Figure 15: Figure S5A: correlation of module eigengenes
ggplot(module_ME_df) +
geom_rect(xmin = -Inf,xmax = Inf,
ymin = -Inf,ymax = Inf,fill="gray90", alpha = 0.1) +
geom_bar(aes(x=id, y=value, fill=sex), stat="identity") +
facet_grid(modules~age, scales="free") +
# theme(axis.text.x = element_text(angle=45, hjust=1)) +
theme(axis.text.x = element_blank(),
axis.title.x = element_blank()) +
scale_fill_manual(values = c("tomato", "steelblue"))
Figure 16: Figure S5B: module eigengenes
module_gene_enrich <- lapply(module_list_genename, function(x) gprofiler(x, organism = "mmusculus", custom_bg = all_detected_genes, hier_filtering = "moderate", min_set_size = 5, max_set_size = 5000, min_isect_size = 3, src_filter = c("GO")))
p1 <- plot_cluster_enrichment(enrich_list = module_gene_enrich[1:4], "CEMi_modules_enrichment_", width=13, exclude_domain = "hp", save = F)
p2 <- plot_cluster_enrichment(enrich_list = module_gene_enrich[5:9], "CEMi_modules_enrichment_", width=13, exclude_domain = "hp", save = F)
p1+p2 +
plot_layout(guides = "collect")
Figure 17: Figure S6C: pathway enrichment result of co-expression modules
save.image("mouse_pituitary_mRNA_analysis.RData")
sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.3 LTS
##
## Matrix products: default
## BLAS: /opt/R/3.6.0/lib/R/lib/libRblas.so
## LAPACK: /opt/R/3.6.0/lib/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_CA.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_CA.UTF-8 LC_COLLATE=en_CA.UTF-8
## [5] LC_MONETARY=en_CA.UTF-8 LC_MESSAGES=en_CA.UTF-8
## [7] LC_PAPER=en_CA.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] pheatmap_1.0.12 CEMiTool_1.8.2
## [3] UpSetR_1.4.0 patchwork_1.0.0
## [5] scatterpie_0.1.4 ggraph_1.0.2
## [7] tidygraph_1.1.2 gProfileR_0.6.7
## [9] fgsea_1.10.0 Rcpp_1.0.1
## [11] dendextend_1.12.0 ggdendro_0.1-20
## [13] cowplot_0.9.4 gridExtra_2.3
## [15] reshape2_1.4.3 RUVSeq_1.18.0
## [17] EDASeq_2.18.0 ShortRead_1.42.0
## [19] GenomicAlignments_1.20.1 SummarizedExperiment_1.14.0
## [21] DelayedArray_0.10.0 matrixStats_0.54.0
## [23] Rsamtools_2.0.0 Biostrings_2.52.0
## [25] XVector_0.24.0 BiocParallel_1.18.0
## [27] ggrepel_0.8.1 ggpubr_0.2.1
## [29] magrittr_1.5 forcats_0.4.0
## [31] stringr_1.4.0 dplyr_0.8.3
## [33] purrr_0.3.3 readr_1.3.1
## [35] tidyr_1.0.0 tibble_2.1.3
## [37] ggplot2_3.2.1 tidyverse_1.2.1
## [39] pcaMethods_1.76.0 Biobase_2.44.0
## [41] gplots_3.0.1.1 RColorBrewer_1.1-2
## [43] rtracklayer_1.44.4 GenomicRanges_1.36.1
## [45] GenomeInfoDb_1.20.0 IRanges_2.18.1
## [47] S4Vectors_0.22.0 BiocGenerics_0.30.0
## [49] edgeR_3.26.5 limma_3.40.2
##
## loaded via a namespace (and not attached):
## [1] ggthemes_4.2.0 R.methodsS3_1.7.1 coda_0.19-2
## [4] ggpmisc_0.3.1 acepack_1.4.1 ffbase_0.12.7
## [7] bit64_0.9-7 knitr_1.33 aroma.light_3.14.0
## [10] R.utils_2.9.0 data.table_1.12.2 rpart_4.1-15
## [13] hwriter_1.3.2 RCurl_1.95-4.12 doParallel_1.0.14
## [16] generics_0.0.2 GenomicFeatures_1.36.4 preprocessCore_1.46.0
## [19] RSQLite_2.1.2 europepmc_0.3 bit_1.1-14
## [22] enrichplot_1.4.0 xml2_1.2.2 lubridate_1.7.4
## [25] assertthat_0.2.1 viridis_0.5.1 xfun_0.29
## [28] hms_0.5.2 jquerylib_0.1.4 evaluate_0.14
## [31] DEoptimR_1.0-8 fansi_0.4.0 progress_1.2.2
## [34] caTools_1.17.1.2 readxl_1.3.1 igraph_1.2.4.1
## [37] DBI_1.0.0 geneplotter_1.62.0 htmlwidgets_1.3
## [40] ellipsis_0.2.0 backports_1.1.4 bookdown_0.11
## [43] annotate_1.62.0 biomaRt_2.40.0 vctrs_0.2.0
## [46] withr_2.1.2 ggforce_0.2.2 triebeard_0.3.0
## [49] robustbase_0.93-5 checkmate_1.9.3 sna_2.4
## [52] prettyunits_1.0.2 cluster_2.0.8 DOSE_3.10.1
## [55] lazyeval_0.2.2 crayon_1.3.4 genefilter_1.66.0
## [58] pkgconfig_2.0.2 labeling_0.3 tweenr_1.0.1
## [61] nlme_3.1-139 nnet_7.3-12 rlang_0.4.2
## [64] lifecycle_0.1.0 modelr_0.1.5 cellranger_1.1.0
## [67] polyclip_1.10-0 graph_1.62.0 Matrix_1.2-17
## [70] urltools_1.7.3 base64enc_0.1-3 ggridges_0.5.1
## [73] viridisLite_0.3.0 bitops_1.0-6 R.oo_1.22.0
## [76] KernSmooth_2.23-15 blob_1.2.0 qvalue_2.16.0
## [79] robust_0.4-18 gridGraphics_0.4-1 ggsignif_0.5.0
## [82] scales_1.0.0 memoise_1.1.0 plyr_1.8.4
## [85] gdata_2.18.0 zlibbioc_1.30.0 compiler_3.6.0
## [88] intergraph_2.0-2 rrcov_1.4-7 cli_2.0.0
## [91] htmlTable_1.13.1 Formula_1.2-3 MASS_7.3-51.4
## [94] WGCNA_1.68 tidyselect_0.2.5 stringi_1.4.3
## [97] highr_0.8 yaml_2.2.0 GOSemSim_2.10.0
## [100] locfit_1.5-9.1 latticeExtra_0.6-28 GeneOverlap_1.20.0
## [103] grid_3.6.0 fastmatch_1.1-0 tools_3.6.0
## [106] rstudioapi_0.10 foreach_1.4.4 foreign_0.8-71
## [109] farver_1.1.0 digest_0.6.19 rvcheck_0.1.3
## [112] pracma_2.2.5 ff_2.2-14 broom_0.5.3
## [115] gRbase_1.8-3.4 httr_1.4.1 AnnotationDbi_1.46.0
## [118] colorspace_1.4-1 rvest_0.3.5 XML_3.98-1.20
## [121] splines_3.6.0 RBGL_1.60.0 ggplotify_0.0.3
## [124] fit.models_0.5-14 xtable_1.8-4 jsonlite_1.6
## [127] dynamicTreeCut_1.63-1 zeallot_0.1.0 R6_2.4.0
## [130] Hmisc_4.2-0 pillar_1.4.3 htmltools_0.3.6
## [133] glue_1.3.1 clusterProfiler_3.12.0 DT_0.7
## [136] DESeq_1.36.0 codetools_0.2-16 pcaPP_1.9-73
## [139] mvtnorm_1.0-11 lattice_0.20-38 network_1.15
## [142] gtools_3.8.1 GO.db_3.8.2 survival_2.44-1.1
## [145] rmarkdown_2.10.6 statnet.common_4.3.0 munsell_0.5.0
## [148] DO.db_2.9 fastcluster_1.1.25 GenomeInfoDbData_1.2.1
## [151] iterators_1.0.10 impute_1.58.0 haven_2.2.0
## [154] gtable_0.3.0